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Abstract 

We address the investigation of the solvation properties of the minimal orientational model for 
water, originally proposed by Bell and Lavis. The model presents two liquid phases separated by 
a critical line. The difference between the two phases is the presence of structure in the liquid of 
lower density, described through orientational order of particles. We have considered the effect of 
small inert solute on the solvent thermodynamic phases. Solute stabilizes the structure of solvent, 
by the organization of solvent particles around solute particles, at low temperatures. Thus, even at 
very high densities, the solution presents clusters of structured water particles surrounding solute 
inert particles, in a region in which pure solvent would be free of structure. Solute intercalates 
with solvent, a feature which has been suggested by experimental and atomistic simulation data. 
Examination of solute solubility has yielded a minimum in that property, which may be associated 
with the minimum found for noble gases. We have obtained a line of minimum solubility (TmS) 
across the phase diagram, accompanying the line of maximum in density (TMD). This coincidence 
is easily explained for non-interacting solute and it is in agreement with earlier results in the 
literature. We give a simple argument which suggests that interacting solute would dislocate TmS 
to higher temperatures. 

PACS numbers: 61.20.Gy,65.20.+w 
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I. INTRODUCTION 



Biological molecules are functional only if organized spatially in very specific arrange- 
ments. This is the case for phospholipids in membranes, proteins soluble in water or mem- 
brane proteins, cholesterol or lipoproteins. One of the main ingredients behind spatial 
organization is solubility: globular proteins maintain their polar moieties on the exterior, in 
contact with water, while membrane protein must turn their polar parts inwards, avoiding 
contact with the hydrophobic bilayer core. 

Solubility depends on chemical structure, but varies with temperature. For simple sub- 
stances, the behavior of solubility with temperature is dependent on miscibility, which de- 
scribes the relative affinities of the molecules in solution The reasoning is simple. If we 
consider the solution phase in equilibrium with the gas phase, two situations exist. Consider 
X to be the solute in solvent Y. If Y and X 'prefer' mixing, which means that the energy of 
an YX pair is lower than the average energy of YY and XX pairs, for the solution energy 
to increase as temperature goes up, X must necessarily leave the solution, thus making sol- 
ubility decrease. On the contrary, if Y and X prefer to phase separate, at low temperatures 
X will go preferentially to the gas phase. However, as temperature goes up, the solution 
energy increases while X dissolves in Y, making the solubility go up. 

Solubility in water is different. Noble gases, for instance, present a temperature of mini- 
mum solubility in water at atmospheric pressure Water presents in numerous thermody- 
namic and dynamic anomalies, and the minimum in solubility is one of them. The origin of 
the anomalies has been investigated theoretically both for statistical and atomistic models. 
However, a simple complete picture has not yet emerged. 

The presence of a hydrogen-bond network was suggested in the '1930s by Bernal js], in 
order to explain the large mobility of H+ and OH- ions: the latter could only be explained if 
protons would jump between neighboring properly oriented molecules in liquid water. The 
idea of an extensive H bond network and a corresponding water structure was probed with X- 
rays for many years, and the presence of the network was confirmed by more recent neutron 
scattering experiments, which pointed to an even more stable structure than previously 
believed 14 1. Hydrogen bonds are considered a key feature in biochemistry [sj. 

The presence of an H-bond network could qualitatively explain the well-known maximum 
in density. The disordering of bonds allows density to increase with temperature since 
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the entropy of the bonds increases while translational entropy decreases, maintaining the 
necessary positive entropy balance. The two entropic effects compete up to a temperature 
at which translational entropy wins over orientational entropy, taking density down, as in 
more 'usual' substances. 

The dynamically connected molecules would also be able to explain the minimum in sol- 
ubility. The contraction of the solvent, driven by decreasing orientational entropy, excludes 
the solute. Thus a decreasing solubility is a consequence of an increasing density of the 
solvent. In this case, the energy of the interactions enters only either to favor the decreasing 
solubility, either to compete with it. 

The study of statistical models capable of displaying properties typical of water has 
led, in the last years, to two basic models: (i) orientational models |6l-ll3|. which reflect 
the H-bonding property of water, and (ii) two-scale isotropic models, inspired on the low- 
temperature low-density property of water. Both models present several of the anomalous 



features of water 
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]2l\ . However, the second kind of model does not involve specific 
orientation of low energy pairs of particles: pair energy is controlled by distance, not by 
orientation. This poses a question of the relevance of the microscopic bonding in relation to 
the macroscopic properties. 

In this study we propose to contribute to further investigation of the relation between 



the solvent structure and solubility. The role of cavity formation in 



;he explanation of 



30 . They examine 



30| . A thorough 



hydrophobic interactions has been recognized by Pratt and Chandler [28 
the difference between cavity formation in associating and simple liquids 
investigation of noble gas solubility was undertaken by Guillot and Guissani 3l| from the 
point of view of atomistic models. Our approach is that of a minimal statistical model. 

shown by us [25|, 127[ to exhibit many anomalous properties in spite of the absence of liquid 
polymorphism. In this study we add non interacting solute particles which occupy a single 
lattice site in order to investigate the effect of solute on solvent properties as well as solute 
solubility. 

The remaining of the paper goes as follows. In sec. II the model without and with solute 
is introduced, simulation details are presented in sec. Ill, the phase diagram of the system 
with solute is shown in sec. IV, the solubility is analyzed in sec. V and conclusions are given 
in sec. VI. 
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II. THE BELL-LAVIS MODEL AS SOLVENT 



The Bell-Lavis (BL) model is defined on a triangular lattice whe^e each site may be empty 
(cTj = 0) or occupied by an anisotropic water molecule (cXi = 1) [7]. Each particle has two 
orientational states, that may be described in terms of six 'arm' variables r*-', with r*-' = 1 
for the bonding state and r*"' = for the inert arm state as illustrated in Fig. [D^a). A 
pair of adjacent molecules interacts via van der Waals with energy e„^^, as well as through 
'hydrogen bonds' of energy ehb, whenever bonding arms point to each other {tI-'t^^ = 1). 
The model is defined by the following effective Hamiltonian in the grand-canonical ensemble 

<i,j> i 

where e^dw and etb are the van der Waals and hydrogen bond interaction energies, respec- 
tively,and ^ is the chemical potential. 

The model phase diagram features depend on the ratio ( = e^dw/^hb (see insets of Fig. [2] 
that illustrates the reduced chemical potential versus reduced temperature for two cases 
of bond strength: weaker, Q = 1/4, and stronger, ( = 1/10). For ( < 1/3, besides the 
gas phase, the model exhibits two liquid phases with different structure. At T = 0, co- 
existence between a gas and a structured liquid of low density (SL) as well as coexistence 
between the structured low density liquid and the non-structured high-density-liquid (NSL) 
are present |32j. However, for finite temperatures, the transition between the two liquids 



becomes critical, as shown from detailed systematic analysis of simulational data 



25|. The 



two liquid phases do not coexist and the density varies continuously at the phase transition 



as shown by susceptibility measurements on sub-lattices density fluctuations [27|. In order 
to stress the absence of a density gap we denominate the two liquid phases as structured (SL) 
and non-structured liquid (NSL), instead of adopting the usual LDL and HDL nomenclature. 
The difference between the two liquid phases lies in the orientational and translational order 
of the bonding particles. The SL phase presents a large population of particles in two of the 
three sub-lattices (see Fig. HJ^b)) associated to a large bonding network, whereas in the NSL 
the density is close to 1 and orientational order is lost. The increase in the temperature 
and the increase in the chemical potential favor the NSL phase. In the case of the stronger 
hydrogen bonds {( = 1/10) the SL is favored and the transition occurs for higher chemical 
potentials. 
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A line of temperatures of maximum density (TMD) lies near the critical line separating 
the two liquid phases. Its pressure and temperature location is not very sensitive to the 
strength of the hydrogen bonds, in both in the ( = 1/10 and in the ( = 1/A cases at low 
chemical potentials is located in the SL phase while for high chemical potentials is located 
at the critical line . In this work we have added inert apolar solutes to the BL model. The 
new particles occupy empty sites and thus interact only via excluded volume with the BL 
solvent particles. Our purpose is the investigation of the effect of the apolar solute upon 
the TMD and the regions of stability of the low and high density phases. Here we address 
the following questions. What would be the effect of adding solute to the structured liquid? 
Under what circumstances does phase separation occur? Is there a solubility minimum? In 
the latter case, can we establish a relation between the density and the solubility anomalies? 



FIG. 1: (a) The arm variables in the triangular lattice, (b) The three sub-lattices, A, B and C. 

III. MONTE CARLO SIMULATIONS 

We have investigated the properties of our model solution through Monte Carlo simula- 
tions, in a mixed ensemble, of fixed chemical potential for solvent and constant density for 
the apolar solute, under periodic boundary conditions. 

The model solvent microscopic configurations were generated through randomly selected 
exclusion, insertion or rotation of water particles, whereas solute movements were based on 
solvent-solute and hole-solute exchanges. Acceptance rates are those of the usual Metropolis 
algorithm: transitions between two configurations are accepted according to the Metropolis 
prescription min{l, exp(— ^A'H)}, where AH is the effective energy difference between the 
two states. Our simulations were carried out for lattice sizes ranging from L = 30 to 
L — 60. Results shown here are for L — 30. All the thermodynamic quantities are expressed 
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in reduced units of ehb and lattice distance. 

IV. SOLVENT PHASE DIAGRAM IN THE PRESENCE OF INERT SOLUTE 

We have investigated how the chemical potential versus temperature phase diagram 
changes by the addition of an inert solute. We study this employing two solute concen- 
trations, 2% and 10%. 

The reduced chemical potential versus reduced temperature phase diagrams for both 
weak and strong bonds, ( = 1/A and C = 1/10, and concentration of 2% of solute, are shown 
in Fig. [21 At this small concentration of solute, the phase diagram suffers small quantitative 
changes: the structures phase SL extends to slightly higher chemical potential, while, in the 
case of weaker bonds, ( = 1/4, the TMD line moves into the SL phase at low temperature. 
Thus solute stabilizes the SL to higher chemical potential, which suggests a reinforcement 
of hydrogen bonding. It also brings down the temperature of maximum density in the case 
of weaker bonds, turning TMD behavior similar to the case of stronger bonds - again, solute 
seems to "strengthen" bonds. 




0.1 0.2 0.3 0.4 0.5 0.1 0.2 0.3 0.4 O.S 0.6 



T T 

FIG. 2: Solvent chemical potential ft vs. reduced temperature T phase diagram for model solution 
at solute concentration 2% and C = 1/4 (left) and C = 1/10 (right). Solid, dashed and dotted lines 
correspond to first-order, second-order phase transitions and the TMD, respectively. The symbols 
C, CE and t denote critical, critical-ending and tricritical points, respectively. The inset displays 
the corresponding phase diagram for pure solvent, for comparison. 

Fig. E] illustrates features of the solution structure in the case of strong bonds, ( = 1/10 
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and concentration of 2% of solute. The triangular lattice is subdivided into three sub-lattice 
as illustrated in Fig. [Hb). The orientation and density of solvent particles, as well as, the 
density of solute are computed on each sublattice. 

The first set of data at Fig. [3]^a)-(c) illustrate the orientation of the solvent molecules, 
the density of solvent and the density of the solute versus temperature for fl = —1.6. The 
graphs show the transition between the NSL to the SL phase by decreasing the temperature. 
It can be seen that in the SL phase (/i = —1.6) solvent occupies mainly two of the sub-lattice 
(with Psoivent nearly 1 at T = 0.3), with complementary orientations (m = -|-1 and m = — 1), 
indicating strong bonding. Both quantities vary abruptly at the transition to the NSL phase, 
around T = 0.45, with homogeneous occupation and orientation of molecules on the three 
sub-lattice. As for solute, at the lower temperature, T = 0.3, occupation of the empty 
sublattice is preferential {psoiute ~ 2%), while the other sub-lattice are nearly empty. As 
solvent disorders on sub-lattice, around T = 0.45, solute densities vary continuously towards 
homogeneous occupations of the three sub-lattice. 

Fig. [3]^e)-(f) illustrate the same data as before but for reduced chemical potential p, = 
—0.4. In this case, no transition is observed. The system is in the NSL phase even at 
low temperature, and sub-lattice solvent orientation and solute density vary continuously 
towards homogeneous distribution on sub-lattice. Solvent density still carries the signature 
of the ordered phase, transitioning smoothly to disorder in a sigmoidal fashion. 

Figure H] displays the reduced chemical potential versus reduced temperature phase dia- 
grams for 10% concentration of the solute, for both values of hydrogen bond strength ( = 1/4 
(left) and ( = 1/10 (right). In this case substantial change in the phase diagrams can be 
seen. The low temperature SL to NSL phase transition seen as one increases chemical po- 
tential for pure solvent is destroyed by the presence of solute. Instead, the transition may be 
reached only from temperature variations, and the SL phase extends to very high chemical 
potentials. The TMD line moves nearer to the critical SL-NSL line and crosses into the SL 
phase. 

In Fig. Owe investigate the solvent orientation and density, as well as, the density of solute 
in each sublattice in different regions of the phase diagram at 10% solute concentration for 
bond strength ( = 1/10. For both low, ft = —1.3, and high, // = 4.0, reduced chemical 
potential, solute orders together with solvent at low temperatures. As can be seen from 
the color identification of lattices, solute goes into the empty lattice while solvent particles 
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FIG. 3: Sub-lattices distributions for solvent and solute properties for 2% solute concentration and 
C = 1/10. From right to left, we have solvent particle orientation mi , solvent particle density pi 
and solute density psi versus reduced temperature T . Top graphs are for lower reduced chemical 
potential p = —1.6, bottom graphs are for higher reduced chemical potential p = —0.4. Colors 
identify sub-lattices. 

orient properly on two sub-lattice in order to connect through bonds. 

For both reduced chemical potentials an abrupt variation of the solvent density, of the 
solute density and of the solvent orientation occur simultaneously, near T = 0.45 for /2 = 
— 1.3 and near T = 0.3 for fx = 4.0 at the SL-NSL transition line. 

However, despite of the qualitative similar behavior there are quantitative important 
differences between the two regions of chemical potential. For the lower chemical potential, 
/i = —1.3, solvent behavior is similar to that of pure solvent. The orientationally ordered 
solvent particles occupy mainly two of the sub-lattice, while the third sublattice remains 
nearly free of solvent. At T = 0.3 nearly 60% of that sub-lattice stands vacant, while solute 
particles occupy 30% of the sites, leaving the other two sub-lattice free of solute. 

At high chemical potential, /i = 4.0, while solute maintains the 30% occupancy of one 
of the sub-lattice at low temperature, the solvent particles fill up the rest of the sublattice 
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FIG. 4: Solvent reduced chemical potential p, versus reduced temperature T phase diagram for 
model solution at solute concentration 10% for C = 1/4 and C = 1/10. Solid, dashed and dotted 
lines correspond to first-order, second-order phase transitions and the TMD, respectively. The 
symbols CE, C and t denote critical-ending, critical and tricritical points, respectively. The insets 
displays the corresponding phase diagrams for pure solvents for comparison. 

sites, reaching 70% occupancy of that sublattice. 

The new behavior induced by the presence of solute is better understood by comparing 
Fig. |3]^e) and Fig. Et^e). Differently from the 2% solute concentration, for the 10% concen- 
tration case the filling up of the lattice yields only partial rupture of hydrogen-bonding. 
Maintenance of the hydrogen bond network at such high density seems to be a result of the 
structuring effect of solute. 

Inspection of typical configurations in different regions of the phase diagrams are quite 
useful at this point. Fig. [7] and Fig. [6] display snapshots of the model system at different 
points (indicated by letters in Fig. H]) in the reduced chemical potential versus reduced 
temperature phase diagrams with solute concentration 10%. 

For C = 1/4 (Fig. ED, at /i = -0.50 and f = 0.25 (inside SL phase, point A in Fig. H), 
lattice is filled up, but patches of structured liquid can be seen with solute localizing only 
in sites which contribute to organize the hydrogen bond network. As the SL-NSL line is 
crossed and for T = 0.50 (point D in Fig. H]), a few isolated solute particles are surrounded 
by water particle structure, while most solute particles are clustered in vacant regions. 

For C = 1/10 (Fig. ED, at /i = -1.40 and f = 0.30 < Ttmd (point C in Fig. HD, a 
fully bonded network of solvent particles is accompanied by solute particles located in the 
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FIG. 5: Solution structure in the new phase. Sub-Lattice solvent particle orientation rm, solvent 
density pi, and solute density psi vs. temperature for different chemical potentials for 10% of solute 
concentration and (" = 1/10. The top and bottom graphs correspond to solvent chemical potentials 
/i = -1.30 and p = 4.0. 



empty sublattice. This gives rise to apparently linear aggregates intercalated by solvent. At 
a temperature higher than the TMD, T = 0.50 (point D in Fig. H]), some bonding of the 
solvent particles in hexagons are still seen, with intercalated solute. However, the system is 
much less dense, and solute particles also localize in large vacant regions. 






FIG. 6: Snapshots of model system for concentration of 10% and C = 1/4 at high chemical potential, 
p, = —0.5, in A (T = 0.25, left) and B (T = 0.5, right) points in Fig. HI Black and Grey circles are 
two orientations of solvent particles, black dots represent solute particles. 
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FIG. 7: Snapshots of model system for concentration of 10% and C = 1 /lO at low chemical potential, 
p, = —1.4, in C (T = 0.3, left) and D (T = 0.5, right) points in Fig. [H Black and gray circles are 
two orientations of solvent particles, black dots represent solute particles. 

V. MODEL SOLUBILITY 



The Ostwald solubiUty S is defined as the ratio between solute densities px in the two 
coexisting phases: 

E = ^. (2) 

Px 

The two coexisting phases, / and //, might either be a gas phase // that coexists with a 
homogeneous hquid phase / [3l| or two Uquid phases / and II, of different relative densities, 
the first poor in solute X, the other reach in X 33 1. 

In the case of liquid-liquid phase separation, the form of the temperature-density co- 
existence curve, at fixed pressure, is indicative of solubility behavior. If the density gap 
decreases as temperature is increased, solubility increases with temperature. However, for 
reentrant coexistence curves, for which the density gap increases as temperature is increased, 
solubility decreases as temperature is raised. 

On the other hand, minimal statistical models show that for dense lattice gas solutions 
with isotropic van der Waals-like interactions solubility behaves univocally with tempera- 
ture. Coexistence densities between an ideal gas mixture and a lattice dense solution, for 
substances Y and X, are obtained from the equality of the corresponding chemical poten- 
tials. Consider interaction constants wyy, wxx and wyx between pairs YY, XX and YX. 
For solute X in the gas phase given by the dimensionless solute X density p^x^ ■, we might 
write 

..gas _ u rpi gas /r,\ 
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whereas for solute X in the dense lattice solution, we have 

^solution ^ ]^ _]l{i-x^f ^ksTln 0;^°'"*^°", (4) 

where w = wyy + wxx — '^wyx and xx is the solution concentration given in mole fraction. 
Equating Eq. © to Eq. (g]) yields 

^solution a 

gas = e 2 e 2 (1 -X ) . (5 
Px 

A slightly different definition of solubility, proportional to the inverse of Henry's constant, 
is given by 



Px IPx 

where p\ is gas density for pure liquid X. Comparing Eq. ([5]) with Eq. ([6]) gives 

= e 2 ) . (7) 



solution 



Thus for poorly miscible solutions, with w < 0, which phase separate at low temperatures, 
solubility increases with temperature, since dY^ jdT oc —w. On the other hand, if the two 
liquids are miscible, when w > 0, solubility decreases as temperature is raised. In either 
case, solubility displays monotonic behavior with temperature. 

The solubility behavior of the dense lattice model is the result of a competition between 
entropy of mixture and an isotropic interaction potential, relies on positional entropy The 
model misses the role of density, an essential feature of water. 

How does the introduction of asymmetry in the interaction potential, accompanied by 
orientational entropy, change this picture? In order to answer to this question we have mea- 
sured the solubility of our model inert solute as a function of temperature for different fixed 
chemical potential of solvent, by assuming coexistence of a gas (phase F) and a homogeneous 
solution phase (X). The gas phase was assumed ideal, thus 

= -fcBTlnpr- (8) 
For the solution phase, the chemical potential of solute was calculated from simulation data 



through Widom's insertion method 3J|. In our semi-grand canonical ensemble the semi- 



grand potential = i^iT,]^, Nx, fJ^) depends on T, V, Nx and on the solvent chemical 

' dip ' 
^dNx' 



potential fi. In the thermodynamic limit, the solute chemical potential = _(-^t) 
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can be approximated by the difference ip^T, V, Nx + I, fi) — ipiT,]^, Nx, /i) and we have 



^s^lution ^ Nx+l] 



(9) 



EiT,V,fi,Nx 

which relates average values in two different ensembles of Nx and Nx+i particles. However, 
the numerator can be interpreted in terms of an average in the ensemble of Nx solute 
particles. Thus we have 

1 \ 



solution 



-ksT In 



-/3A« 



f^pSolution\^ 



(10) 



where Am is the additional energy due to insertion of solute molecule to a system of solute 



concentration px 
solubility 



olution 



Nx/V. Finally, by equating and /i^*, we obtain for the 



^solution 

Px 



-l3Au 



V,Nx,l^- 



(11) 



In the Fig. [8] we display our data for solubility S versus temperature for bond strength 
( = 1/10. As can be seen, a minimum is present for different chemical potential of solvent. 
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FIG. 8: Ostwald coefficient S versus T for different /j, for ( = 1/10 and 10%. The black line 
corresponds the temperature for which the density presents a maximum. 



The temperature of minimum solubility (TmS) coincides entirely with the temperature 
of maximum density (TMD), in the /2 vs. T plane. This is to be expected for inert solutes. 
Inspection of Eq. (fTTj) for inert solutes yields exp{— /JA-u} = 1 for insertion into empty 
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sites and 0, otherwise. Thus solubihty can be directly related to the overall liquid density 

{Nx + N'"^''^'^^)/V. Thus 

E/„—0Au\ -I solvent solute /lO\ 

= (e ^ )T,V,Nx,f^ = -L - P - Px 5 U^j 

and for fixed solute density 

dT ~ ^^^^ 
and therefore the TMD is accompanied by the TmS. Is the coincidence between TmS and 
TMD restricted to inert solutes? 

It is tempting to extend our analysis of Eq. f lTT]) to interacting solutes. A first simplest 
approach to the question would be to investigate the energetic effect on solubility through 
the following approximation 

V — ~ n ^solvent solute\ -l3<Au> (TA^ 

thus 



dT 



dp'"'"''''' ^ _ ^solvent _ psolute^ ( < > _ < Au > 



dT ' rA / , ^^^2 - ^J' 



e-^<'^">. (15) 



Since < Am > is necessarily negative and '^'^^^ is necessarily positive, this result implies 
that the minimum in solubility should occur at a temperature higher than TMD. This is in 
accordance with data on solubility of gases. 



VI. CONCLUSIONS 



In this paper we have considered the investigation of the thermodynamic phases and of 
solubility of the Bell-Lavis (BL) water model in the presence of small inert solute. The 
Bell-Lavis two-dimensional orientational model presents a density anomaly and two liquid 
phases of different structure^. 

We have considered two fixed concentrations of solute, respectively 2% and 10%. In both 
cases, but more evidently for 10% solute, the presence of solute 'strengthens' the hydrogen 
bonds. Inspection and comparison of phase diagrams show that the structured phase is 
stabilized to higher temperatures, at fixed chemical potential, and to much higher chemical 

^ a modified form of the BL model fias been investigated as to solvation entropy and enthalpy properties [ssl 

15 



potentials. For the higher concentration of solute, the transition between the structured and 
the unstructured phases as chemical potential is varied disappears. Examination of solvent 
structure shows that the presence of solute nucleates patches of hydrogen-bonded solvent 
particles. Solute intercalates with properly oriented solvent. Both features, increment of 
water structure and solvent-separated solute states, have been reported from experiments 
and atomistic models |36l |. 



Solubility of our small inert solutes presents a minimum (TmS), which coincides with 



the maximum solvent density (TMD), as expected [28|. For interacting solute, a simple 
argument leads us to expect TmS to occur at higher temperatures for the BL solvent model. 

Investigation of the latter point, as well as the effect of solute size are the subject of 
ongoing work. 
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